#include "ewpt.h"
#include <cstring>
#include <iostream>
#include <sstream>
#include <cmath>
#include <vector>
#include "TMath.h"

//defined as mass*sqrt(2)/v
double EWPT::top_couple=EWPT::top_mass*sqrt(2)/v;
double EWPT::W_couple=EWPT::W_mass*sqrt(2)/v;
double EWPT::Z_couple=EWPT::Z_mass*sqrt(2)/v;
double EWPT::h_couple=EWPT::h_mass*sqrt(2)/v;

//constructor
EWPT::EWPT()
{
  //initialize other functions
  V_tree=new TF1("V_tree", "[0]*(x^4-2*x^2*[1])", -5000, 5000);

  double quartic=h_mass*h_mass/(2*v*v);
  double v_phi=v/sqrt(2.);
  double v_phi2=v*v/2;

  V_tree->SetParameter(0, quartic/2.);
  V_tree->SetParameter(1, v_phi2);

  //string for mass factors
  //g*v^2
  string mass_v="([2]*[1])";
  string mass_h="([2]*x^2)";

  //coleman weinberg
  double factor=1/(64*PI*PI);

  stringstream formula;
  formula << "[0]*("
    
    //log term 2*mh^4*log(mh/mv)
	  << "2*"<<mass_h<<"^2*log("
	  <<mass_h<<"/"<<mass_v<<")"

    /*
    //-2*(mh^2-m_v^2)(3mh^2-m_v^2)
	  << "-2*(" <<mass_h<<"-"<<mass_v<<")"
	  << "*(3*" <<mass_h<<"-"<<mass_v<<"))";
    */

    //-3/2 mh^4
	  << "-3*"<<mass_h<<"^2/2"

    //+2 mh^2 mv^2
	  << "+2*"<<mass_h<<"*"<<mass_v<<")";
    

  V_loop_top=new TF1("V_loop_top", 
		     formula.str().c_str(),
		     -5000, 5000);

  //dirac fermion = 4 degrees of freedom
  V_loop_top->SetParameter(0, -4*factor);
  V_loop_top->SetParameter(1, v_phi2);

  //top yukawa coupling
  V_loop_top->SetParameter(2, top_mass*top_mass/v_phi2);


  //now for gauge bosons
  V_loop_W=new TF1("V_loop_W", 
		   formula.str().c_str(),
		   -5000, 5000);
  V_loop_W->SetParameter(0, 2*3*factor);
  V_loop_W->SetParameter(1, v_phi2);
  V_loop_W->SetParameter(2, W_mass*W_mass/v_phi2);

  //for Z boson
  V_loop_Z=new TF1("V_loop_Z", 
		   formula.str().c_str(),
		   -5000, 5000);
  V_loop_Z->SetParameter(0, 3*factor);
  V_loop_Z->SetParameter(1, v_phi2);
  V_loop_Z->SetParameter(2, Z_mass*Z_mass/v_phi2);
  
  //for higgs
  V_loop_H=new TF1("V_loop_H", 
		   formula.str().c_str(),
		   -5000, 5000);
  V_loop_H->SetParameter(0, factor);
  V_loop_H->SetParameter(1, v_phi2);
  V_loop_H->SetParameter(2, h_mass*h_mass/v_phi2);
  
  //test
  double x=200;
  /*
  cout<<"Tree: "<<V_tree->Eval(x)<<endl;
  cout<<"top: "<<V_loop_top->Eval(x)<<endl;
  cout<<"W: "<<V_loop_W->Eval(x)<<endl;
  cout<<"Z: "<<V_loop_Z->Eval(x)<<endl;
  cout<<"h: "<<V_loop_H->Eval(x)<<endl;
  */

}

double EWPT::SM(double x, double T) const
{
  //first part = zero temperature
  double V0= V_tree->Eval(x);

  if(x*x < 1e-20)
    x=1e-10;
  
  V0+=
    V_loop_top->Eval(x) + V_loop_W->Eval(x) +
    V_loop_Z->Eval(x) + V_loop_H->Eval(x);


  //now temperature dependent terms
  if(T>0)
    {
      //we want to compute mass/T for different particles
      //first T^4
      double T4=T*T*T*T;

      //m over T neglecting coupling
      double mT=x/T;

      //multiplicity factors:
      //top = 4, W=6, Z=3, h=1
      V0+=T4*(
	      4*tint_fermion(top_couple*mT) +
	      6*tint_boson(W_couple*mT) +
	      3*tint_boson(Z_couple*mT) +
	      1*tint_boson(h_couple*mT)
	      );
    }
  
  /*
  cout<<"Tree: "<<V_tree->Eval(x)<<endl;
  cout<<"top: "<<V_loop_top->Eval(x)<<endl;
  cout<<"W: "<<V_loop_W->Eval(x)<<endl;
  cout<<"Z: "<<V_loop_Z->Eval(x)<<endl;
  cout<<"h: "<<V_loop_H->Eval(x)<<endl;
  */

  return V0;
}

TGraph* EWPT::GetGraph(double T) const
{
  vector<double> x;
  vector<double> y;
  for(int i=-300; i<=300; i+=2)
    {
      x.push_back(i);
      y.push_back((*this)(i,T));
    }

  TGraph* g=new TGraph(x.size(), &x[0], &y[0]);

  stringstream sout;
  sout<<"T_"<<T;
  g->SetName(sout.str().c_str());
  return g;
}

double EWPT::smooth(double r) const
{
  if(r<0) return 1.0;

  else if(r>1) return 0.0;

  return r*r*(3.-2.*r);
}

double EWPT::tint_boson(double r) const
{
  if(r>100)
    return 0;

  const double factor=-1.0/(2*PI*PI);

  //make sure r>0
  r=fabs(r);

  //if r>matching_high, use the Bessel function approximation
  if(r>matching_high)
    return factor*
      (bessel_sum(r)+bessel_sum(2*r)/16+bessel_sum(3*r)/81)
      - b0;

  double r2=r*r;
  double r3=r2*r;
  double r4=r2*r2;
  double r4log=0;
  if (r>0) r4log=r4*log(r);
  double r6=r4*r2;

  //if r < matching_low, use high temp expansion
  if(r<matching_low)
    return b2*r2 + b3*r3 + b4*r4 + b4log*r4log + b6*r6;

  //if r between matching_low and matching_high, use the transition region
  
  double p_smooth= 
    smooth((r-matching_low)/(matching_high-matching_low));

  return 
    p_smooth*
    (factor*
     (bessel_sum(r)+bessel_sum(2*r)/16+bessel_sum(3*r)/81)
     - b0) +
    (1-p_smooth)*
    (b2*r2 + b3*r3 + b4*r4 + b4log*r4log + b6*r6);

}


double EWPT::tint_fermion(double r) const
{
  if(r>100)
    return 0;

  r=fabs(r);

  const double factor=-1.0/(2*PI*PI);

  //if r>matching_high, use the Bessel function approximation
  //subtract the zero point energy
  if(r>matching_high)
    return factor*
      (bessel_sum(r)-bessel_sum(2*r)/16+bessel_sum(3*r)/81)
      - f0;

  double r2=r*r;
  double r4=r2*r2;
  double r4log=0;
  if (r>0) r4log=r4*log(r);
  double r6=r4*r2;

  //if r < matching_low, use high temp expansion
  if(r<matching_low)
    return f2*r2 + f4*r4 + f4log*r4log + f6*r6;

  //if r between matching_low and matching_high, use the transition region
  
  double p_smooth= 
    smooth((r-matching_low)/(matching_high-matching_low));

  return 
    p_smooth*
    (factor*
     (bessel_sum(r)-bessel_sum(2*r)/16+bessel_sum(3*r)/81)
     -f0) +
    (1-p_smooth)*
    (f2*r2 + f4*r4 + f4log*r4log + f6*r6);

}


TGraph* EWPT::GetIntGraph(int spin) const
{
  vector<double> x;
  vector<double> y;

  for(double i=0; i<=20; i+=0.2)
    {
      //cout<<"calculating : "<<i<<endl;
      x.push_back(i);
      
      if(spin==1)
	y.push_back(tint_fermion(i));

      else
	y.push_back(tint_boson(i));
    }

  TGraph* g=new TGraph(x.size(), &x[0], &y[0]);
  g->SetName("int");
  return g;
}


TGraph* EWPT::GetLowGraph(int spin) const
{
  vector<double> x;
  vector<double> y;

  for(double i=0; i<=20; i+=0.2)
    {
      //cout<<"calculating : "<<i<<endl;
      x.push_back(i);

      double& r=i;
      double r2=r*r;
      double r3=r2*r;
      double r4=r2*r2;
      double r4log=0;
      if (r>0) r4log=r4*log(r);
      double r6=r4*r2;

      
      if(spin==1)
	y.push_back(f0 + f2*r2 + f4*r4 + f4log*r4log + f6*r6);

      else
	y.push_back
	  (b0 + b2*r2 + b3*r3 + b4*r4 + b4log*r4log + b6*r6);
    }

  TGraph* g=new TGraph(x.size(), &x[0], &y[0]);
  g->SetName("int_low");
  return g;
}

TGraph* EWPT::GetHighGraph(int spin) const
{
  vector<double> x;
  vector<double> y;

  const double factor=-1.0/(2*PI*PI);

  for(double i=0; i<=20; i+=0.2)
    {
      //cout<<"calculating : "<<i<<endl;
      x.push_back(i);

      double& r=i;
      
      if(spin==1)
	y.push_back(factor*
	(bessel_sum(r)-bessel_sum(2*r)/16+bessel_sum(3*r)/81));

      else
	y.push_back
	  (factor*
	   (bessel_sum(r)+bessel_sum(2*r)/16+bessel_sum(3*r)/81)
	   );
    }

  TGraph* g=new TGraph(x.size(), &x[0], &y[0]);
  g->SetName("int_high");
  return g;
}

double EWPT::Dx(double x, double T, double step) const
{
  return ((*this)(x+step, T) - (*this)(x-step, T))/(2*step);
}

double EWPT::Dx2(double x, double T, double step) const
{
  return ((*this)(x+step, T) + (*this)(x-step, T) - 
	  2*(*this)(x, T))/(step*step);
}

double EWPT::find_v(double T, double step, double guess) const
{
  //try for 100 times maximally
  //start at 500
  double x=guess;
  bool success=false;

  for(int i=0; i<100; ++i)
    {
      double inc= Dx(x, T, step)/Dx2(x, T, step);
      if(fabs(inc) > 1e4)
	{
	  cout<<"ERROR: infinite step size, no vev found"<<endl;
	  x=0;
	  break;
	}

      x-= inc;

      if(fabs(inc)<2*step)
	{
	  success=true;
	  break;
	}
    }

  if(success)
    return x;
  

  cout<<"ERROR: minimum not found after 100 iterations, returning zero"<<endl;
  return 0.0;

}


double EWPT::Tc(double& vev) const
{
  //cout<<"Tc called"<<endl;
  //first some initialization
  double T=0;
  vev=find_v(0);
  
  //maximum 100 times iteration
  bool success=false;
  //initial step size
  double step=10.;
  for(int i=0; i<100; ++i)
    {
      //cout<<"iteration: "<<i<<endl;

      double oldT=T;
      double oldvev=vev;
      Tc_helper(T, vev, step);
      
      if(fabs(T-oldT) < 1e-3 && vev>1e-3 &&
	 fabs(vev-oldvev) < 1e-3 && step < 1e-4)
	{
	  success=true;
	  break;
	}
      step/=3.;
    }
  
  if(!success)
    {
      cout<<"ERROR: cannot find Tc, return 0"<<endl;
      vev=0;
      return 0;
    }

  return T;
}

void EWPT::Tc_helper(double& T, double& vev, double step) const
{
  //cout<<"initial T,vev: "<<T<<", "<<vev<<endl;

  if(vev<1e-3) return;
  //now start with temperature T at vev
  //find the T that still has a nonzero vev
  for(double Tp=T+step; Tp<10000; Tp+=step)
    {
      double vevp=find_v(Tp);
      if(vevp < 1e-3 || (*this)(vevp, Tp)>0 )
	{
	  // cout<<"final T,vev: "<<T<<", "<<vev<<endl;
	  return;
	}

      //update the vev if there is a higher temperature
      //with a non-zero vev and the field value is lower
      vev=vevp;
      T=Tp;
    }
}
